﻿#include <learn/math/integrator.h>

using namespace xi;
using namespace xi::math;

Integrator::Integrator(std::size_t N) : m_N(N)
{
}

std::size_t Integrator::node_num() const
{
    return m_N;
}

Integrator_Trapezoidal::Integrator_Trapezoidal(std::size_t N) : Integrator(N)
{
}

double Integrator_Trapezoidal::operator()(double (*f)(double), double a, double b)
{
    std::size_t N = m_N - 1;

    double h = (b - a) / N;
    double sum = (f(a) + f(b)) / 2;
    for (std::size_t i = 1; i < N; i++)
        sum += f(a + i * h);

    return sum * h;
}

Integrator_Simpson::Integrator_Simpson(std::size_t N) : Integrator(N)
{
}

double Integrator_Simpson::operator()(double (*f)(double), double a, double b)
{
    std::size_t N = m_N - 1;
    N = N / 2 * 2;

    double h = (b - a) / N;
    double sum = 0.0;
    for (std::size_t k = 1; k < N / 2 + 1; k++)
    {
        double r = a + (2 * k - 2) * h;
        sum += f(r) + 4 * f(r + h) + f(r + 2 * h);
    }

    return sum * h / 3;
}

Integrator_Gauss::Integrator_Gauss(std::size_t N) : Integrator(N)
{
    MatrixXd Tn(N, N);
    Tn.set_zero();
    for (std::size_t i = 1; i < N; i++)
        Tn(i, i - 1) = Tn(i - 1, i) = 1 / std::sqrt(4 - 1.0 / (i * i));

    Matrix_SymEigen_Solver solver(std::move(Tn));
    solver.perform();
    auto V = solver.solve(N * 10);
    auto D = solver.move_matrix();

    for (std::size_t i = 0; i < N; i++)
    {
        nodes.push_back(D(i, i));
        weights.push_back(2 * V(0, i) * V(0, i));
    }
}

double Integrator_Gauss::operator()(double (*f)(double), double a, double b)
{
    int N = weights.size() - 1;

    double sep = (b - a) / 2;
    double mid = (a + b) / 2;
    std::vector<double> y(N + 1);
    for (std::size_t i = 0; i < N + 1; i++)
        y[i] = sep * f(sep * nodes[i] + mid);

    return std::inner_product(weights.begin(), weights.end(), y.begin(), 0.0);
}